Anthropogenic events and responses to environmental stress are shaping the genomes of Ethiopian indigenous goats

Anthropological and biophysical processes have shaped livestock genomes over Millenia and can explain their current geographic distribution and genetic divergence. We analyzed 57 Ethiopian indigenous domestic goat genomes alongside 67 equivalents of east, west, and north-west African, European, South Asian, Middle East, and wild Bezoar goats. Cluster, ADMIXTURE (K = 4) and phylogenetic analysis revealed four genetic groups comprising African, European, South Asian, and wild Bezoar goats. The Middle Eastern goats had an admixed genome of these four genetic groups. At K = 5, the West African Dwarf and Moroccan goats were separated from East African goats demonstrating a likely historical legacy of goat arrival and dispersal into Africa via the coastal Mediterranean Sea and the Horn of Africa. FST, XP-EHH, and Hp analysis revealed signatures of selection in Ethiopian goats overlaying genes for thermo-sensitivity, oxidative stress response, high-altitude hypoxic adaptation, reproductive fitness, pathogen defence, immunity, pigmentation, DNA repair, modulation of renal function and integrated fluid and electrolyte homeostasis. Notable examples include TRPV1 (a nociception gene); PTPMT1 (a critical hypoxia survival gene); RETREG (a regulator of reticulophagy during starvation), and WNK4 (a molecular switch for osmoregulation). These results suggest that human-mediated translocations and adaptation to contrasting environments are shaping indigenous African goat genomes.


Genome-wide genetic diversity and population structure
The mean values of the genetic diversity indices are shown in Table 1.Among Ethiopian goats, the mean and standard deviation of H O , H E , π, and D ST were 0.347, 0.371, 0.0021, and 0.303, respectively.The H O , H E, and D ST were highest in Ethiopian goats.The F HOM and F RoH averaged 0.064 ± 0.04 and 0.059 ± 0.01, respectively in Ethiopian goats; the highest values were in Italian (F HOM = 0.246 ± 0.16) and Pakistan (F HOM = 0.243 ± 0.11) goats.In Ethiopian goats, the mean number and length of RoH were 932.9 ± 90.31 and 174.14 ± 36.31Mb, respectively; the highest values were in Pakistan goats and averaged 1663.4 ± 811.17 and 410.27 ± 233.37 Mb, respectively.
Genetic structure and relationship were assessed at the first instance for Ethiopian indigenous goats only and then for the overall dataset that included all Ethiopian, non-Ethiopian and the wild Bezoar goats.PC1 and PC2 of the PCA explained 3.67% and 2.78% of the total genetic variation of Ethiopian goats (Fig. 1a) and each revealed two broad genetic clusters.The first cluster identified by PC1 includes Afar, Hararghe Highland, Shorteared Somali, Long-eared Somali, and Woyto-Guji.The second cluster of PC1 comprises Abergelle, Gonder, Agew, Ambo, Gumuz, Arsi-Bale and Keffa.Similarly, PC2 grouped together Afar, Hararghe highland, Shorteared Somali, Abergelle, Gonder, Agew, Ambo, Gumuz, and Arsi-Bale in its first cluster while its second cluster grouped together Long-eared Somali, Woyto-Guji, and Keffa.In combination, these two PC's reveal four clusters that likely reflect a fine-scale genetic divergence in Ethiopian indigenous goats.The first one is made up of Keffa; the second comprises Long-eared Somali and Woyto-Guji; the third comprises Afar, Hararghe Highland and Short-eared Somali, and the fourth includes Abergelle, Gonder, Agew, Ambo, Gumuz and Arsi-Bale.
To further examine the genetic structure of Ethiopian indigenous goats, we generated an NJ phylogeny using F ST genetic distances (Fig. 1d).This revealed two broad genetic clusters that support the PCA results.We named these two genetic clusters A and B. The F ST phylogeny also supports the fine-scale and deep genetic structure in Ethiopian goats revealed by PCA with Keffa being genetically distinct.ADMIXTURE analysis showed the lowest CV error at K = 1 suggesting genetic homogeneity among Ethiopian goats.However, the genetic profile at K = 2 (Fig. 1b) shows two broad genetic clusters, whose population composition corresponds to the ones identified by PC1 of the PCA and the F ST phylogeny.For brevity, we also name these A and B. At the same K value (i.e., K = 2), ADMIXTURE analysis provides additional insights into the genetic structure of Ethiopian goats that are not apparent in the PCA and F ST phylogeny; that the genomes of Arsi-Bale, Woyto-Guji and Hararghe Highland comprise different proportions of clusters A and B. Cluster A predominates in Arsi-Bale (69.0%) and, B in Woyto-Guji (71.47%) and Hararghe Highland (80.47%) (Supplementary Fig. S5).An analysis of the geographic

Genome-wide dynamics
Genome-wide demographic dynamics were inferred by assessing LD patterns against genomic distances and changes in N e over generation time for each population and genetic clusters that were revealed by PCA and ADMIXTURE.The pattern of LD decay was the same for all the population clusters (Figs.1e, 2e).It reveals higher LD at shorter genomic distances which decays rapidly reaching a plateau around ~ 0.2 Mb.Generally, Ethiopian goats show higher average LD (r 2 ≥ 0.25) (Fig. 1e) and within them, Gonder had the highest LD (r 2 ≥ 0.5), followed by Ambo and Abergelle (r 2 ≥ 0.34) and then the other populations (r 2 ≥ 0.25).
The pattern of N e was similar across all the Ethiopian indigenous goat populations (Fig. 1f).There was a gradual decline in N e between 1000 and 500 generations ago, following which the decline accelerated.Performing the analysis based on the genetic groups generated by cluster analysis reveals that both clusters A and B of  Ethiopian populations show an increase in N e between 1000 and 400 generations ago following which there is a drastic decline (Fig. 2f).

Selection signature analysis
We investigated genome-wide selection signatures with Hp, F ST, and XP-EHH tests to explore whether the population clusters observed in Ethiopian goats are the result of adaptation to different environments.For this analysis, we selected Afar, Arsi-Bale, and Keffa goats as they occurred in different clusters on the PCA.The H P test revealed a total of 196 candidate regions under selection, that overlapped 484 genes across Afar, Arsi-Bale, and Keffa goats (Fig. 3a; Supplementary Table S5).The F ST (Fig. 3b; Supplementary Table S6) and XP-EHH (Fig. 3c; Supplementary Table S7) tests identified a total of 222 and 356 candidate regions, respectively that spanned 411 and 757 genes when contrasting Afar and Arsi-Bale, Afar and Keffa, and Arsi-Bale and Keffa.Based on the ARS1 RefSeq gene annotation, a total of 145 regions did not overlap with any gene(s) and/or spanned genes that are not yet annotated (Supplementary Tables S5-S7).Given the large number of candidate regions identified by the three tests, to retain high specificity we used the threshold score values of − 7.0 (ZH P ), 7.0 (ZF ST ), and 6.0 (XP-EHH) to define the top-most significant selection regions.We regarded these to be the primary selection signatures that are shaping the genomes of the study populations.Our results and discussions will be based on these top-most significant regions unless specified otherwise.Based on the above cut-off threshold scores, eight candidate regions in Afar (overlapping 36 genes), seven in Arsi-Bale (seven genes), and eight in Keffa (55 genes) were identified by ZH P as the top-most significant regions (Table 2).The ZF ST identified nine regions in Afar vs Arsi-Bale (13 genes), nine in Afar vs Keffa (34 genes), and four in Arsi-Bale vs Keffa (4 genes) (Table 3) while those identified by XP-EHH were, 23 for Afar vs Arsi-Bale (80 genes), seven for Afar vs Keffa (12 genes), and nine for Arsi-Bale vs Keffa (12 genes) (Table 4).Of these 253 genes, eight (SCNNIB, SMPD3, NUP43, ZNF609, COG7, ZFP90, PCMT1, PIFI) overlapped between ZH P and ZF ST, one each (Fig. 4) overlapped between ZH P and XP-EHH (BPIFB4), and ZF ST and XP-EHH (ALOX5AP).
Functional annotation and gene ontology analysis were performed with DAVID at two levels (1) for all the genes identified by ZH P in each population and (2) for all the genes identified by the three comparative analyses (Afar vs Arsi-Bale, Afar vs Keffa, and Arsi-Bale vs Keffa) involving ZF ST and XP-EHH combined as one gene list.We found several functional gene clusters that were significantly enriched (Supplementary Tables S8 and  S9).In Afar, the enriched clusters were RNA degradation, endocytosis, autophagy, Rap1 signaling pathway, PI3K-Akt signaling pathway, and MAPK signaling pathways.In Arsi-Bale, enrichments were Ras signaling pathway, Rap1 signaling pathway, PI3K-Akt signaling pathway, MAPK signaling pathway, bacterial invasion of epithelial cells and melanogenesis.In Keffa, MAPK signaling pathway, glycerophospholipid metabolism,TGFbeta signaling pathway, Yersinia infection and IL-17 signaling pathway were overrepresented.In the Afar vs Arsi-Bale comparative analysis, significantly enriched clusters were for neurological system processes, apoptotic processes involved in morphogenesis, and sensory organ development.In Afar vs Keffa, the significantly enriched

Discussion
Following goat domestication from the wild Bezoar ~ 11,000 years ago 1,2 , population bottlenecks, inbreeding, intermixing, and selection (natural and artificial) have been modifying the genomes of domestic goats.This is especially the case in Africa where their initial introduction as exotics from the centre of domestication into the continent exposed them to novel environments and physiological extremes and challenges.Here, by analysing genetic divergence by exploiting SNPs from 57 individual genomes from 12 Ethiopian indigenous domestic goat populations and placing them within a comparative dataset of published African and non-African domestic goat breeds and the wild Bezoar goat allowed us to contextualize regional, continental, and intercontinental diversity and divergence.Notwithstanding differences in sequencing depth and platforms, our sequence statistics (Supplementary Table S1) were consistent with previous observations in goats [18][19][20] , sheep 21,22 , and cattle 23,24 , indicating the high quality and reliability of our dataset.
Irrespective of the species and markers used, several studies have reported high genetic diversity estimated from whole-genome sequences in indigenous livestock compared to exotic/commercial breeds 22,25,26 .The reduced genetic diversity in the latter is the outcome of long-term artificial selection and/or genetic drift due to a demographic history of low effective population sizes.The genome-wide autosomal SNPs showed high values in all the estimated parameters of genetic diversity in Ethiopian goats (Table 3).These values are comparable with those reported for indigenous goats in Uganda 27 , South Africa 28 , Egypt 14 , Spain 29 and Italy 30 estimated using SNP microarray genotypes and from whole-genome sequence analysis in indigenous African goats 25 , cattle 24 and sheep 22 .In diploid genomes, RoH represents continuous homozygous segments of DNA sequences and can provide insights into how population history, structure, and demography have evolved over time 31 .In this study, estimates of F RoH and F HOM reveal low levels of inbreeding in all the populations analysed which is consistent with findings in Ugandan 32 and Egyptian Barki goats 14 .Our analysis also revealed a similar pattern of high frequency (> 50%) of the shortest average RoH length segments (0.1-0.25 Mb length category) in the populations analysed (Supplementary Table S4; Supplementary Figs.S3, S4).This skewed distribution agrees with other findings on goats 32,33 , sheep 34 and cattle [34][35][36] where long RoH segments were more infrequent than shorter ones.Short ROHs suggest inbreeding is not recent.Taken together, these results appear to suggest that the high heterogeneity in the study populations could be the result of a combination of factors.These include low inbreeding due to random mating and historic admixture arising from the communal use of resources and/or the www.nature.com/scientificreports/common practice of sharing and/or gifting stock to cement social bonds and relationships, and the predominance of natural selection which favours standing genetic variation.This is supported by the demographic dynamics, which show all populations have low LD ranges (> 200 kb) and historic high N e (Fig. 2).The highest average LD was recorded in Ethiopian (0.25-0.55) goat populations.These LD values are within the range reported for a large number of goat breeds 37 and East African shorthorn Zebu cattle 38 .However, when the data is combined the lowest and highest LD values are recorded in Ethiopian and wild Bezoar goats, respectively.The low levels of long-range LD show the lack of intensive selection or the populations have had large effective ancestral population sizes 39 .The decline in N e from around 1250 years ago (500 generations) is difficult to explain as no significant events that could have affected goat populations have been reported in the Horn of Africa region.We however speculate that the start of this decline could have been driven by the Lapanarat-Mahlatule drought that occurred in the region and was characterized by a sequel of severe droughts and political upheavals 40 .
Though not native to Africa, goats are ubiquitous in the continent and are closely associated with the subsistence practices, socio-cultural life, and economy of many African societies.Archaeological evidence suggests that goats first arrived in Africa from Southwest Asia 41 .Our phylogenetic and population structure analyses at the whole genome level revealed at least two ancestral genetic groups in African goats; one was observed in west and north-west African goats and the other in east African goats.We hypothesise that these two genetic ancestries could be contiguous with the trajectories of the initial movement of the first goat-farming pastoralists into the continent from Southwest Asia.Radiocarbon dates of caprovid remains from the North African Mediterranean coastline are amongst the oldest in the continent 6,42 .Thus, the west/north-west African ancestry, due to its location furthest from the postulated initial entry point of goats into Africa, was likely the first in the continent.To arrive at their current locations this ancestral group may have dispersed along two routes following its entry via Egypt.The ancestry that occurs in north-west African goats dispersed along Africa side of the Mediterranean Sea while the one found in west Africa dispersed overland across the present-day Sudano-Sahelian belt 43 .
The East African genetic ancestry could represent goats that spread from the Near East into the central Sahara, Sudan Nile, and the Ethiopian highlands between 6500 and 5000 BP 6,7 .This genome ancestry could have spread  1a,b,d), and their respective geographic distribution across the country could shed light on the dispersal of this genomic ancestry into East Africa.The two clusters mirror findings from the analysis of 50 K SNP genotype data of Ethiopian goats 15 and of mitochondrial DNA, which identified two haplogroups in the Ethiopian 11 , Kenyan 10 , Sudanese 44 and Egyptian 45 indigenous goats.Cluster A predominates in goats found in the north and west of Ethiopia, while B occurs at a higher frequency in goats found in the south and east of the country (Fig. 1d).Their spread across the country was most likely facilitated by socio-cultural and commercial interactions as can be inferred from anthropologic, linguistic and human genetic studies 46 .This geographic spread led us to suggest that cluster A could have dispersed to Ethiopia from Egypt following the Nile River basin or across the Red Sea Hills, while cluster B most likely arrived in Ethiopia via the Horn of Africa through the Bab el-Mandeb strait.Worth mentioning is that the inclusion of genomes from east, west, and north-west Africa, Europe, South Asia, the Middle East, and from the wild Bezoar goats provided an interesting insight.The ADMIXTURE analysis showed that all the four genomes it revealed are present in Iranian (Middle East) goats, which supports this as the cradle of present-day goat genome diversity.
Ethiopia is characterized by a diverse combination of agro-eco-climates and ancient and modern human ethnic diversity that may have influenced the genome architecture of indigenous livestock.Our phylogenetic analysis revealed fine-scale genetic structuring in Ethiopian goats, which we hypothesize could be driven by environmental adaptation.Thus, signatures of selection were investigated within and between Afar, Arsi-Bale, and Keffa on the premise that their divergence is driven by adaptation to contrasting environments and therefore can serve as good proxies for investigating selection signatures resulting from genomic divergence.Afar goats inhabit a low altitude area (120-200 masl) with a hot semi-arid/arid agro-ecology (mean annual rainfall 150-300 mm).Arsi-Bale goats inhabit the Bale mountains (> 3000 masl) with a cool, cold sub-humid and alpine agro-ecology.These two populations clustered separately in the PCA suggesting genomic divergence.Keffa, which also occurs in a mid-altitude (≤ 1800 masl) environment, showed a clear genetic divergence from the other Ethiopian populations, suggesting it is genetically unique.This was observed previously by Tarekegn 11 from the analysis of 50 K SNP genotype data.The selection signature analysis revealed several candidate genomic regions under selection, some of which did not span any genes.This is not uncommon; it has been reported in cattle 47 , sheep 14 , and goats 15 .In line with our hypothesis, the top-most candidate regions under selection, spanned genes with roles in adaptation to different biophysical stressors rather than production.This suggests that natural rather than deliberate artificial selection is the principal driver of divergence in the studied populations and is mediated by www.nature.com/scientificreports/ the concurrent action of a complex network of genes.The large number of candidate regions and genes detected is also not surprising.Similar results have been reported for livestock species from extreme environments 21,48 .
Our findings corroborate those of Tarekegn et al. 15 who reported the divergence of Keffa from other Ethiopian goats by analysing 50 K SNP genotype data suggesting it to be genetically distinct.The authors speculated that it was due to Keffa being trypanotolerant.Whereas, we found a large number of selection signals in Keffa, Tarekegn et al. found none.We attribute this difference to the higher resolution afforded by whole-genome sequences.The XP-EHH between Afar vs Keffa detected a strong selection region on CHI9 (73.95-74.85Mb) spanning two genes, PPP1R14C and IYD.Interestingly, PPP1R14C was found in a selection sweep region in the Ethiopian Sheko cattle and it was postulated to be one of the candidate genes contributing to trypanotolerance in the breed 49 .Both Sheko cattle and Keffa goats occur in an area where trypanosomosis is an economically important livestock disease.Whether this signature presents convergent evolution in the two species for trypanotolerance is difficult to say from the current data.Cattle and goats are both bovids with minor genomic differences 50,51 indicating minor genomic differences since their divergence from a common ancestor 14 .Therefore, it is not unusual to find genetic and biological similarities between the two species.While this signature could contribute to the uniqueness of Keffa goats, we cannot conclude it is the only factor.Afar goats reside in a semi-arid and arid environment.It is characterized by complex interacting biophysical stressors including heat, physical exhaustion, direct solar radiation, and resource (feed and water) scarcity.It is thus unsurprising that the selection tests for Afar goats, revealed signatures that spanned genes for thermosensitivity e.g., TRPV1, PCMT1, PACRGL, and IYD.The activity of PCMT1 was reported to peak under lethal temperatures, suggesting a role for the enzyme in short-term responses to heat extremes 52 .D' Alessandro et al. 53 demonstrated that PCMT1 also plays an essential role that ensures normal RBC circulation during oxidative stress.TRPV1, a thermally activated ion channel, plays a major role in thermosensation, thermoregulation, and nociception 54 , and its activation provides a high-temperature noxious-heat-avoidance signal 55,56 .We also found strong selection signatures in regions overlaying ANGPTL7 and ENKD1, which are associated with maintaining skin integrity.By regulating extracellular matrix formation coupled with its high expression in keratinocytes 57 , ANGPTL7 plays a central role in skin homeostasis, repair, and regeneration 58 .ENKD1 regulates spindle orientation in basal keratinocytes, which promotes epidermal stratification and regeneration 59 .The epidermis protects an organism from extremes of the external environment, reduces water and heat loss and pathogen entry 60 .Thus, in arid environments, skin integrity is important for protection against solar radiation, maintaining internal homeostasis and moisture loss.
In several of the top-most candidate regions in Afar, we also found five genes, FGF23, TIGAR , RAMP2, SLC12A3, and WNK4, which are critical in regulating mineral and nutritional metabolism and homeostasis, and water balance.FGF23 regulates mineral homeostasis and plays an essential physiological role in phosphate and vitamin D metabolism 61 .Bensaad et al. 62 showed that TIGAR can lower intracellular reactive oxygen species in response to nutrient starvation or metabolic stress, and functions to inhibit autophagy.The adrenomedullin (AM) and its receptor-modulating protein, RAMP2, facilitate early adaptation to cardiovascular stress by maintaining and regulating cardiac mitochondria and cardiovascular homeostasis against cardiovascular stress 63 .This is critical in countering physical stress resulting from long-distance trekking in search of pasture and water which can put pressure on cardiomyocytes.The AM-RAMP2 system also suppresses endoplasmic reticulum stress-induced tubule cell death, thereby exerting a protective effect on kidneys 64 .Kidneys are critical in water and electrolyte (Na + , Cl − , and K + ) homeostasis.SLC12A3 encodes the thiazide-sensitive sodium chloride cotransporter, which is primarily expressed in the kidney, intestines, and bones.It plays a key role in sodium, potassium, and blood pressure regulation in response to various hormonal and non-hormonal stimuli 65 .WNK4 functions as a molecular switch that varies the balance between NaCl reabsorption and K + secretion to maintain integrated homeostasis 66 .It plays a key role in coordinating the activities of flux pathways, which are regulated by the renin-angiotensin-aldosterone system, to achieve integrated fluid and electrolyte homeostasis (osmoregulation) in the distal nephron and distal colon 67,68 following acute dehydration and rapid rehydration.
Arsi-Bale and Keffa goats are found in Bale and Keffa zones, respectively.Keffa goats live between 1000 and 1800 masl, while Arsi-Bale goats inhabit elevations between 3000 and 4377 masl.These zones are characterized by high levels of precipitation.High-altitude environments impose a selective constraint in the form of hypobaric hypoxia 69 , which results in insufficient oxygen supply in body tissues.This would affect normal physiological www.nature.com/scientificreports/functions and can result in organ failure and death 70,71 .Our selection signature analyses for Arsi-Bale and Keffa goats revealed several strong signatures that spanned genes such as RUNDC3B, TIGAR , PTPMT1, STXBP4, and ALOX5AP relating to hypoxia adaptation.It has been shown that RUNDC3B is amongst genes with variants associated with increased risk of high-altitude polycythemia (HAPC) in Tibetan dwellers 72 .HAPC is characterized by excessive proliferation of circulating erythrocytes due to high-altitude hypobaric hypoxia.The compensation for prolonged hypobaric hypoxia exposure is the main reason for the change in erythrocyte production and haemoglobin concentration that elevates oxygen retention, transportation, and exchange.Findings by Kimata et al. 73 showed that TIGAR is a significant mediator of cellular energy homeostasis (glycolysis) and cell death (apoptosis) under ischemic/hypoxic stress.Through a genome-wide CRISPR-Cas9 KO library screening, Bao et al. 74 identified PTPMT1, an important enzyme for cardiolipin synthesis, as the third most significant gene for hypoxic adaptation/survival, ranking right after HIF-1α and HIF-1β.In a genome-wide study of genetic adaptation to high altitude in feral Andean Horses, the highest region of allele frequency divergence spanned, amongst five other genes, STXBP4 and COX11 75 .A highly significant association was also found between these two genes and are in strong LD in both humans and horses.COX11 is up-regulated in chronic hypoxia suggesting a role in dealing with oxygen deficit, possibly by acting as a heme biosynthetic enzyme that transports copper to heme 76,77 .Its strong association and LD with STXBP4 may suggest a similar function for this gene.ALOX5AP was identified in sheep as a potential candidate for climate-mediated adaptation 78 .In humans, a mutation in ALOX5AP was associated with lung function 79 .Given the high altitude and restricted oxygen concentration in the Ethiopian highlands, ALOX5AP may play a role in adaptation by modulating respiratory function.
Extreme environments (high-altitude, semi-arid, and arid) impose anatomical, physiological, and metabolic challenges with strong evolutionary pressure due to long-term exposure to acute and chronic stress and other factors dependent on the natural history of a population.These factors are an underlying constant in the three test populations (Arsi-Bale, Keffa, Afar) and what is likely driving their difference is their exposure to biotic and abiotic stress factors in contrasting environments.It is therefore insightful that some of the strongest selection signals spanned genes for oxidative stress mitigation (TRPM2, DDX28, ALDH3B1, TIGAR , ALOX5AP, RAMP2, POU2F1, HSF4), DNA damage repair and maintenance of genome stability (EXOSC10, CTCF, FZD5, PIF1, FBRSL1, HSF4), skin pigmentation and characteristics (FBRSL1, EZH1, POU2F1, SOX5, KITLG, OAZ2, HSF4) and reproduction function (NUP43, EXOSC10, TARDBP, DPEP3, ESCO2, OAZ2, HSD17B1, PSMC3IP, EZH1, NECTIN3, KATNAL1, CACNB2, ELF5, HASPIN, USP42, SUN5, KITLG, KCTD19, TSNAXIP1).Long-term exposure to acute and/or chronic stress results in an imbalance in the production and accumulation of free radicals (reactive oxygen species) which can induce oxidative stress.Oxidative stress causes base damage and DNA strand breaks resulting in apoptosis and necrosis.Therefore, oxidative stress response and DNA strand repair can preserve genome integrity and normal mechanisms of cellular signaling under stress.
In conclusion, despite the complexity of the agro-ecological and climatic conditions of Africa, domestic goats occur across the continent.The results presented here show at least two genomic ancestries in African goats, and that genetic divergence and genomic plasticity is the driver of the successful integration of the species into African environments.These findings are significant in the context of improving livestock productivity in the continent in view of the projected consequences of climate change on biodiversity.A complete characterization of African indigenous goat's unique genomic variation and adaptation can inform the formulation and design of breeding programs that promote the long-term sustainable goat productivity.

DNA samples and sequencing
Twelve indigenous Ethiopian goat populations (Supplementary Table 1; Fig. 5b,c) that were previously described in Tarekegn 11,15 were used in this study.No ethics permissions were required.DNA samples of five individuals per population were selected at random and whole genome sequenced in Novogene, China (https:// en.novog ene.com/ servi ces/ reser achse rvices/ genome-seque ncing/ whole-genome-seque ncing-wgs/).Whole-genome sequencing was performed on an Illumina NovaSeq 6000 Platform (Illumina, San Diego, CA, USA) and 150 bp of paired-end reads at a target coverage depth of 10× were generated.The quality of the sequences was assessed with FASTQC v0.11.5 80 and three samples failed quality control (Supplementary Table S1).For comparative genome analysis and referencing, whole-genome sequences of 67 individuals of the east (Kenya, Boran goat), west (Nigeria, Dwarf goat) and north-west (Morocco) African (n = 20), South Asian (Pakistan and Bangladesh; n = 17), Middle East (Iran, n = 10) and European (Italy; n = 10) domestic goats, and 10 of wild Bezoar goats obtained from public databases (https:// www.goatg enome.org/ vargo ats.html) were included in the study (Supplementary Table S1; Fig. 5).The wild Bezoar goats, an extant species found in western Asia from Turkey to Pakistan is the presumed ancestor of modern-day domestic goats 2 .All clean reads were aligned to the C. hircus reference genome (ARS1; GenBank accession number GCA_001704415.1) using the BWA tool v0.7.17 81 .The alignment files in SAM format were converted to BAM format with SAMtools 82 .
We applied the GATK v3.8 workflow 83 for variant calling and discovery.SAMTools was used to sort and index the alignment files and duplicate reads were purged using PICARD v2.18.2 (https:// broad insti tute.github.io/ picard/).Base Quality Score Recalibration was used to create recalibrated BAM files from the recalibrated table.Variant Quality Score Recalibration was performed using VariantRecalibrator in GATK with the "knownSites" set to the C. hircus dbSNP reference panel (https:// e99.ensem bl.org/ capra_ hircus).Individual VCF files were called using the HaplotypeCaller in GATK followed by GenotypeGVCFs to generate consolidated GVCF files that contained the raw SNPs and Indels.Variant refinement was performed using variantRecalibrator in GATK.Variant call annotations such as Read Depth, Quality of Depth, Fisher Strand Test, Mapping Quality Score, Mapping Quality Rank Sum Test, Read Position Rank Sum Test Statistic, StrandOddsRatio Test, mode SNP and the VQSRTranchesSNP90 to 100 were used to produce recalibrated SNPs.Using ApplyRecalibration, a tranche sensitivity threshold of 99% was used to generate filtered variants, and post-processing using SelectVariants was conducted to remove variants that failed the GATK filtering parameters.In this step, a total of 26,990,415 and 56,648,064 markers remained which, included autosomal and sex chromosomes for Ethiopian, and other (African, Eurasian, and wild Bezoar) goats, respectively.Further, quality control using the command line 'SelectVariant' and 'restrictAlleles' were used to remove sex chromosomes and multiallelic SNPs, respectively.Finally, the subsequent analyses were performed using 24,759,579 and 42,728,409 autosomal biallelic markers found across 57 and 124 individuals from 12 and 25 populations, respectively (Table 1).

Genome-wide genetic diversity and dynamics
Within-population variation was determined by estimating observed (H O ) and expected (H E ) heterozygosity, nucleotide diversity (π), within-population genetic distance (D ST ), and two inbreeding coefficients (F HOM and F RoH ).The H O , H E and π were estimated with VCFTools v0.1.15 84.To compute π, a 20 kb window size and a sliding step of 10 kb was used.For each population and group of populations, SNP density was also calculated with VCFtools v0.1.15with the command line "-SNPdensity1000" and its mean and standard deviation was computed with R v4.1.0 85.The D ST was calculated using the command line "-genome" in PLINK v1.9 86 .The genetic distance between all individuals within a population was then calculated as D = 1 − D ST .
Demographic dynamics were investigated by assessing pairwise LD between all pairs of autosomal biallelic variants (SNPs) over genomic distance through the correlation coefficient (r 2 ) with PLINK v1.9.Historic

Figure 1 .
Figure 1.Population genetic structure and relationship of the Ethiopian goat populations based on (a) PCA and (b) ADMIXTURE analysis at K = 2, (c) geographic distribution and genetic admixture proportion of the Ethiopian indigenous goat populations (d) phylogenetic tree constructed using F ST values, (e) the pattern of linkage disequilibrium (r 2 ) from 0 to 1 Mb and (f) the pattern of effective population size (N e ) in the past 1000 generations.

Figure 4 .
Figure 4. Venn diagram showing the number of genome-wide selection signals in Ethiopian indigenous goats that are specific to, and shared between H P , F ST, and XP-EHH approaches.

Figure 5 .
Figure 5. Map of the study areas (a) The geographic location of the African, South Asian, Middle Eastern, European, and Bezoar goat populations, (b) the geographic distribution of Ethiopian goat populations based on elevation, and (c) Agro-ecological zones (Sources: Own processed maps using global data set in ArcGIS environment version 10.8).

Table 1 .
Indicators of genetic diversity and variation in the 25 Ethiopian and non-Ethiopian goat populations analysed in the study.Π = Nucleotide diversity, H O = Observed heterozygosity, H E = Expected heterozygosity, D ST = Genetic distance, F HOM and F ROH = Genomic inbreeding coefficients, ROH length = length of Run of homozygosity (RoH), SD = Standard deviation, Mb = Mega base,

Table 2 .
Candidate regions revealed to be under selection and associated genes as identified by the ZHp approach in Afar, Arsi-Bale and Keffa Ethiopian indigenous goat populations.

Table 3 .
Candidate regions under selection and associated genes as identified by the ZF ST approach in Afar x Arsi-Bale, Afar x Keffa and Arsi-Bale x Keffa Ethiopian indigenous goat populations.

Table 4 .
Candidate regions under selection and associated genes as identified by the XP-EHH approach in Afar x Arsi-Bale, Afar x Keffa and Arsi-Bale x Keffa Ethiopian indigenous goat populations.